program main2x2
	
	implicit NONE
	
	integer(4), parameter			:: L = 100
	integer(4), parameter			:: ND = 16
	integer(4)						:: Hzeeman(ND,ND)
	integer(4)						:: Hexchng(ND,ND)
	
	real(8)							:: Hm(ND,ND)
	
	integer(4)						:: hblck1(4,4), hblck2(4,4), hdg(ND)
	
	integer(4)						:: ip, iq
	integer(4)						:: sts
	
	real(8)							:: Jex, hx
	real(8)							:: evals(ND), E(ND), Z(ND), E2(ND)
	real(8)							:: emin
	
!!!!
	integer(4)						:: lwork = 3*ND
	real(8), allocatable			:: wrk(:)
	
	integer(4)						:: eunit
	real(8)							:: t
	integer(4), parameter			:: tmax = 300
	real(8)							:: tt
	real(8)							:: Zs, Es, Esqs, Eavr, Esqavr, Cavr
!!!!
	
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
	open( newunit = eunit, file='evals.txt', status = 'replace', action = 'write')
	
	hblck1 = 0
	forall( ip = 1:4) hblck1(ip,ip) = 1
	hblck2(:,1) = [ 0, 1, 1, 0]
	hblck2(:,2) = [ 1, 0, 0, 1]
	hblck2(:,3) = [ 1, 0, 0, 1]
	hblck2(:,4) = [ 0, 1, 1, 0]
	hdg = [	4, 0, 0, 0, 0,-4, 0, 0, 0, 0,-4, 0, 0, 0, 0, 4]
!!!!
	Hexchng = 0
	forall( ip = 1:ND) Hexchng(ip,ip) = hdg(ip)
!!!
	Hzeeman = 0
	forall( ip = 1:4, iq = 1:4)
		Hzeeman(ip,iq) = hblck2(ip,iq)
		Hzeeman(4+ip,iq) = hblck1(ip,iq)
		Hzeeman(8+ip,iq) = hblck1(ip,iq)
		
		Hzeeman(ip,4+iq) = hblck1(ip,iq)
		Hzeeman(4+ip,4+iq) = hblck2(ip,iq)
		Hzeeman(12+ip,4+iq) = hblck1(ip,iq)
		
		Hzeeman(ip,8+iq) = hblck1(ip,iq)
		Hzeeman(8+ip,8+iq) = hblck2(ip,iq)
		Hzeeman(12+ip,8+iq) = hblck1(ip,iq)
		
		Hzeeman(4+ip,12+iq) = hblck1(ip,iq)
		Hzeeman(8+ip,12+iq) = hblck1(ip,iq)
		Hzeeman(12+ip,12+iq) = hblck2(ip,iq)
	end forall
!!!
	Jex = 1.0d0
	hx = 0.0d0
	Zs = 0.0d0
	Es = 0.0d0
	Esqs = 0.0d0
	
!	do ip = 1, ND
!		print '(20I2)', Hzeeman(:,ip)
!	end do
	
	Hm = - hx * Hzeeman + Jex * Hexchng

	allocate( wrk(lwork), stat = sts)
	call DSYEV( 'N', 'U', ND, Hm, ND, evals, wrk, lwork, sts)
	deallocate( wrk)
	
	if ( jex > 0.0d0) then
		emin = 5 * jex
	else
		emin = - 5 * jex
	end if
	
!	emin = minval(evals)
	evals = evals - emin
	
	do tt = 1, 40
!	do tt = 100, 100
		
		t = tt * 0.1d0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		forall(ip = 1:ND) Z(ip) = expSerL(- evals(ip) / t, L)
		E(:) = evals(:) * Z (:)
		E2(:) = evals(:)**2 * Z (:)
		
		Zs = sum( Z(:))
		Es = sum( E(:))
		Esqs = sum( E2(:))
	
		Eavr = Es / Zs
		Esqavr = Esqs / Zs
		Cavr = ( Esqavr - Eavr**2) / t**2
		
		Eavr = Eavr + emin
	
!		write (eunit, *) t, Eavr, Cavr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		forall(ip = 1:ND) Z(ip) = expSerL(- evals(ip) / t, L)
		Zs = sum( Z(:))
		forall(ip = 1:ND) Z(ip) = expSerL(- evals(ip) / t, L-1)
		E(:) = evals(:) * Z (:)
		E2(:) = evals(:)**2 * Z (:)
		
		Es = sum( E(:))
		Esqs = sum( E2(:))
	
		Eavr = Es / Zs
		Esqavr = Esqs / Zs
		Cavr = ( Esqavr - Eavr**2) / t**2
		
		Eavr = Eavr + emin
	
!		write (eunit, *) t, Eavr, Cavr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
		Z(:) = exp(- evals(:) / t)
		E(:) = evals(:) * Z (:)
		E2(:) = evals(:)**2 * Z (:)
		
		Zs = sum( Z(:))
		Es = sum( E(:))
		Esqs = sum( E2(:))
	
		Eavr = Es / Zs
		Esqavr = Esqs / Zs
		Cavr = ( Esqavr - Eavr**2) / t**2
		
		Eavr = Eavr + emin
	
		write (eunit, *) t, Eavr, Cavr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
	end do
	
	continue
	
	close (eunit)
	
contains
	
pure function expSerL( x, LL)
	! Variables
	real(8)					:: expSerL
	real(8), intent(IN)		:: x
	integer(4), intent(IN)	:: LL
	!----------------------------------------------
	integer(4)			:: ip
	real(8)				:: y, z

	! Body of expSerL
	
	y = 1.0d0
	z = 1.0d0
	do ip = 1, LL
		z = z * x / ip
		y = y + z
	end do
	expSerL = y	

end function expSerL

end program main2x2